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Abstract. This paper discusses the exact simulation of the stock price 
process underlying the 3/2 model. Using a result derived by Craddock 
and Lennox using Lie Symmetry Analysis, we adapt the Broadie-Kaya 
algorithm for the simulation of affine processes to the 3/2 model. We 
also discuss variance reduction techniques and find that conditional 
Monte Carlo techniques combined with quasi-Monte Carlo point sets 
result in significant variance reductions. 
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1 Introduction 



Exact simulation allows us to sample solutions of stochastic differential equations 
(SDEs) from the appropriate distribution functions. Alternatively, one could 
discretize the time interval and simulate the solution by stepping through the 
time grid. The latter approach has two drawbacks. Firstly, a bias is introduced, 
on which it is often difficult to obtain a priori estimates. Secondly, the time 
discretization usually increases the computational complexity of the simulation. 
Regarding the topic of time discretizations for SDEs, we refer the reader to [20] , 
and for methods which show how to trade-off bias and variance, we refer the 
reader to [I2] , and also the recent work on multilevel methods, see for example 
and [13j. Exact simulation methods enjoy the appealing properties of avoiding a 
simulation bias and the increased computational complexity due to time-stepping. 
Finally, they also allow us to asses the quality of discretization schemes. 

For some SDEs, the correct distribution of the solution is well-known or easy 
to obtain, such as Brownian motion and its direct transformations, such as geo- 
metric Brownian motion and the Ornstein-Uhlenbeck process. Furthermore, for 
processes in the family of the squared Bessel process, such as the Bessel and the 
square- root process, exact simulation schemes are also known, see [28] and [27] . 
Recently, the exact simulation of general diffusion processes has become topical, 
see for example the works [2], [3], [1], and [7]. 

In a related paper, [5] solved an important problem of both practical and theoret- 
ical interest, namely the exact simulation of stochastic volatility and affine jump- 
diffusion processes, in particular, the Heston model and its extensions. Roughly 
speaking, the method developed in [5j can be described as follows: One firstly 
samples the integrated variance conditional on which the stock price follows a log- 
normal distribution, whose variance parameter is determined by the integrated 
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variance. To be more precise, [5] first sample the variance at the final time 
point, using the known distribution of the solution of the square-root process, 
and consequently derive the Laplace transform of the conditional distribution of 
the integrated variance, relying on a result from [26] . 

In this paper, we adopt an analogous approach: we also sample the variance at 
the final time point first, and consequently derive the Laplace transform of the 
conditional distribution of the integrated variance. The technique is analogous 
to the one employed by [5] and is found in |9j , where Lie Symmetry Methods are 
employed to derive Laplace transforms of functionals of diffusions, such as squared 
Bessel processes. Having obtained the Laplace transform of the distribution, we 
have reduced the problem to sampling from the lognormal distribution, which is 
trivial. 

The 3/2 model was introduced in [16], and also studied in [23], [6], and [H], and 
is empirically supported, see [18] . It is a stochastic volatility model, and is related 
to the Heston model, [T5|, in the following way: Under the Heston model, the 
variance process is modeled via a square-root process, under the 3/2 model, the 
variance process is modeled via the inverse of a square-root process. We point out 
that the variance process underlying the 3/2 model, the inverse of a square- root 
process, has also been used to model interest rates, see e.g. [1]. 

From a mathematical point of view, the problem is also interesting: though not 
an affine process, the 3/2 model is still analytically tractable, in particular, the 
characteristic function of the logarithm of the stock price still has a closed-form 
solution. Invoking a result obtained via Lie Symmetry Analysis, which also allows 
us to recover the result presented in [26|, we manage to obtain the Laplace trans- 
form of the conditional distribution of the integrated variance. This emphasizes 
that results from Lie Symmetry Analysis can also be employed when designing 
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Monte Carlo methods, whereas so far the main apphcation of Lie Symmetry Anal- 
ysis to finance has been in the derivation of closed-form pricing formulae, see e.g. 

We also discuss variance reduction techniques to speed up the simulation. Since 
[29] . conditional Monte Carlo techniques have been known to improve the effi- 
ciency of Monte Carlo algorithms, especially when combined with quasi-Monte 
Carlo methods, see also [5j. The results we present confirm this observation: we 
manage to achieve considerable variance reductions. 

The remainder of the paper is structured as follows. In Section [21 we state the 
algorithm we employ to simulate the stock price process under the 3/2 model. 
The main challenge of the algorithm lies in the simulation of the conditional 
integrated variance, which is discussed in Section [3l In Section HJ we detail the 
implementation of the algorithm and show some numerical results. These results 
are improved on in Section [5l where variance reduction techniques are studied. 
Section [6] concludes the paper. 

2 The Monte Carlo Algorithm 

The dynamics of the stock price under the 3/2 model under the risk-neutral 
measure are given by 

7 Q 

= rdt + ^/VtpdWl + v^Vl - pHW^ , (2.1) 
dVt = KVt{e -Vt)dt + eVt^^^dWt\ (2.2) 

Equation (12. ip describes the dynamics of the stock price, and Eq. (12. 2p the 
dynamics of the variance process. We denote by and W'^ two independent 
Brownian motions. Regarding the parameters, r represents the constant interest 
rate, p the instantaneous correlation between the return and the variance process 
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and e governs the volatility of volatility. The speed of mean reversion is given 
by K,Vt and 6 denotes the long-run mean of the variance process. However V = 
{Vt , t > 0} is just the inverse of a square-root process, as we now demonstrate, 
in particular, we introduce the square-root process X = {Xt , t > 0} via Xt = ^. 
The dynamics of Xt are given by 

dXt ={K + e^- neXt) dt - e^/xldW^^ . 

Hence, using the process X, we obtain the following dynamics for the stock price, 
where u > t, 

Su = Stexp!^r{u~t)-^j^ {X.y'ds + pj^ (v^j^'rfiy^j 
exp|v/r^^'(v^)"'diy2| 

We aim to adopt the approach from [5]. In this regard, it is useful to study 
log(Xt), for which we obtain the dynamics 

d\og{Xt)= (^^±-^-Kd^ dt~e{yxi^'^ dWl. 

Hence 

log(X„)=log(Xi)+(^«:+0^ ^-K9{u-t)-eJ^ (v^)"'rfH^^ 
or equivalently 

i\^y - 7 ('"^ (^) ^ /" I - '^"c - «)) • p-^) 

We hence arrive at Algorithm^ which is analogous to the Broadie-Kaya algorithm 
from [5]. 
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Algorithm 1 Exact Simulation Algorithm for the 3/2 model 



Step 1) Simulate X„ 



Xt using the noncentral x^-distribution 



Step 2) Simulate /" ;g Xt,Xu 

Step 3) Recover {^/Xy^ dW} from Eq. Q 

Step 4) Simulate Su given St, (V^I) ^ dW} and (A^) ^ ds via 



log(5'„) 



1 



A ( log(5t) + r(M - t) - I (A,)-^ ds + p^ 



where 



0" 



Regarding Step 1) of Algorithm [H from the dynamics of A it is clear, see e. 
[T9] . that conditional on A^, 



A„ exp {k6'('U — t)} 2 



c{u — t) 



where 6 = a = -rzi\i and 



c(t) = (exp {K^t} - 1) /{AkO) . 

Step 2) is discussed in Section [3l in fact, we derive the Laplace transform of the 
conditional distribution of 

Xt, Xu ■ 



t Xs 



Steps 3) and 4) are trivial. 
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3 The Distribution of the Conditional Integrated 
Variance 



We now discuss the simulation of y' 
from 151. Recall that 



Xt.Xu, where we follow the approach 



K + e 



- ds + j (-e) ^sdWl . 



We firstly change the constant in front of y/X, to 2: 



-e) ^sdWl 



'X_. 







{-dWl) 



4 



2 / ^/XJW^2 
'o 

2 



Xu.2.,^dWl2^. 

V \~^J 4 



2 2 

Now, setting u = we get du = ^ds, and hence 



ill 
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^{^) = ^0 + ^ y \^ ~X^\du + 2 

Defining ^{u) = Xau we arrive at 



X^dWl 



We now introduce n = '^'^^^2^"^^ and j = — then 



m = m + f +n)du+2 f ^^i{^)dwl 

Jo Jo 



(3.4) 



To develop a formula for 



E{exp^-a I ^ 





Xj- . X, 



we firstly apply a change of law formula to eliminate the random component of 
the drift term and get a process with j = in (13. 4p . i.e. a squared Bessel process. 
Again, we proceed as in [5]. 

ds 



E exp < —a 



Xt, Xu 





\ 4a f 











where the last equahty follows from formula (6.d) in [26] and E denotes the 
expectation taken with respect to the law of the squared Bessel process. So we 
need formulae for 



E \exp< — - 



Aa f 4, dl f 



m 2 74 



17 /. 



and 



E I exp <( I 



eh 



(3.5) 



(3.6) 



respectively. Regarding Eq. fl3.6p . we have the following result, see e.g. 



E(exp<j-^ /" as)ds 



bt 



sinh{bt) 



exp 



x + y 
2t 



1 - btcoth{bt)) \ 



Iu{b^/xy / sinhibt)) 



where z/ = ^ — 1, the index of the squared Bessel process. Regarding Eq. (13. 6p . 
we proceed as follows: 

"* ds 



E(exp{-aat)-^- j i{s)ds-c 



exp {—ay} E ( exp 





^ Jo 



as) 



^{s)ds — c 



ds 




^{o) = x,m = y 



where ql'^'* {x, y) denotes the transition density of started from a; at being in y 
at time t. Consequently, we can look at 



£^ exp <^ -a^{t) 



u2 rt 



2 Jo 



C,{s)ds — c 



ds 
W) 



as the Laplace transform of 
52 rt 



* ds 



e(0)=x,e(t) = y)gf^(x,i/) 



This Laplace transform is explicitly inverted in [S], see Example 5.3: 



i{(,)^x,i{t)^y\qf\x,y) 



and 



where v = n/2 — 1. Hence we have 

ds 



E exp < —a* 



t 



T ( jVXtXu ' 

\/!^2+8a/(e2)'^ sinh(j A) ' 

T ( jyxtXu \ 

-''^Vsmh(jA)/ 



(3.7) 



where A = ^ — We have characterized the conditional distribution of 



u ds_ 

t Xs 



Xt,Xu- We will sample from this distribution by inversion, which we dis- 



cuss in Section m 



4 Implementation 



In this section, we discuss the implementation of Algorithm [H where we rely on [5j 
and |14j. The implementation of Step 1) was already discussed in Section [2l and 
Steps 3) and 4) are trivial, so we focus on Step 2). We obtain the characteristic 
function of 

Xt , Xu 



by setting a* = —ia in (13. 7p and we use the notation 



$(a) = E (^exp !^ia ^| 



Xt, Xu 
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Furthermore, we define the conditional distribution function 



F{x) := Pi r^<x 



2 1"°° sin{ux) 



Xt,X^] = - / — ^ — '-^{u)du 



TT ./n U 



where we do not emphasize the dependence of the cumulative distribution func- 
tion on Xt and Xu- We employ the trapezoidal rule to approximate the integral 
in Eq. (14.81) numerically, and obtain 

X hx 2 ^ sin(/i?x) „ .,, 

where /i is the grid size associated with the trapezoidal rule and denotes the 
number of terms in the summation; we use ed{h) to denote the discretisation error 
associated with grid size h and cxiN) denotes the truncation error resulting from 
the termination of the sum after N terms. From the discussion in [5J, it is known 
that choosing 



where 1 — F{ue) = e and < x < u^, results in a discretization error ed{h) of 
magnitude e. However, as is difficult to obtain in this way, we set equal to 
the mean plus 12 standard deviations. Regarding the value of A^, as in [5J, we 
terminate the sum at j = N, where 

\^{hN)\ ^ TTt 



N 2 
and e is the desired discretization error. 

Having discussed the implementation of the conditional cumulative probability 
function, F{x), we proceed as follows. Following [14], Section 4.3, we evaluate 
F(-) on the grid 

2 — 1 

Xi = + ^^{u, - wfi),i = 1, . . . ,M + 1 , (4.8) 
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where denotes the conditional expected value of the integrated variance and we 
choose w = 0.01 and M = 200. Consequently, we sample from the conditional 
distribution by inversion. 

Clearly, the most time consuming step when sampling from the conditional dis- 
tribution is the evaluation of the modified Bessel function of the first kind, Iv{z), 
which has to evaluated at complex v. We point out that this operation is eas- 
ily done in MATHEMATICA, and consequently we use the MATHEMATICA 
computing package to perform the simulation. We remark that we also perform 
the inversion of the probability distribution using MATHEMATICA, unlike [5], 
where the equation, 

F{x) = U, 

where U is simulated from a uniform [0, 1] distribution, was solved for x us- 
ing Newton's method. The reason we do not employ their technique is twofold: 
Firstly, [5] could use a good initial guess for Newton's method, as they could ap- 
proximate the conditional distribution of the integrated variance by the Gaussian 
distribution. The analogous result for the 3/2 model is not known, but we found 
that the Newton search is highly dependent on the initial guess. Secondly, using 
MATHEMATICA's built-in functions, we can perform this search both quickly 
and reliably, in particular, we compute the probability distribution on the grid 
of points dSHD, and use MATHEMATICA's NEAREST function to identify the 
point at which the probability distribution assumes a value closest to the simu- 
lated uniform random variable. We point out that the computational time taken 
is highly dependent on one's experience with MATHEMATICA's procedural pro- 
gramming, and should the modified Bessel function of the first kind, allowing 
for V to be complex, become available for other computing packages, or should 
one choose to implement it, the computational times can be expected to differ 
substantially. Hence we report the number of sample paths required to achieve 
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Number of simulation trials 


Price estimates 


Standard error estimates 


2560 


0.46787175 


0.02721607 


10240 


0.430909 


0.0132623 


40960 


0.442416 


0.00672314 



Table 1: Price estimates and standard error estimates for a European call option 
using Monte Carlo simulation 

a particular standard error. This metric is platform and user independent, so 
we find it useful to present. Furthermore, it ties in nicely with the variance re- 
duction techniques presented in Section [5l which allow us to identify the number 
of trajectories by which the computational effort is reduced. Finally, the Monte 
Carlo algorithm we present is of course parallelizable, in principle, one trajectory 
could be computed on one processor, which allows us to substantially reduce the 
computational effort, and we expect this trend to continue in the future. 

We conclude this section by applying Algorithm [T] to the pricing of a European 
call option. This product is chosen, as it allows us to verify Algorithm [T] using a 
different method, namely we price the European call option via the characteristic 
function of the logarithm of St- The characteristic function of the logarithm of 
5*7- is well-known, see e.g. [6], Theorem 3, or [16], [23]. We choose the following 
set of parameters 

So = l, K =1, K = 2,9 = 1.5,e = 0.2, p = -0.5 , T = 1 , r = 0.05 (4.9) 

and obtain the reference value 0.443059. Table [T] shows price estimates and 
standard errors for a given number of simulation trials. 
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5 Variance Reduction Techniques for Stochastic 
Volatihty Models 



We recall some well-known variance reduction techniques for stochastic volatility 
models, see in particular |29]. These methods are not restricted to the 3/2 model, 
but are more generally applicable, see also j5]. The key observation is that given 
paths of 

Vsds and [ ^/VsdW^ , 
Jo 

the price of a European call option is given by the Black-Scholes price, with 
modified initial share price 

5o = 5oexp|-^^ Vsds + pj^ ^/V.dWl 

and adjusted volatility cta/T^-^, where 




Consequently, using BS {Sq, K, r, r, a) to denote the Black-Scholes price of a Eu- 
ropean call with initial stock price Sq, strike K, interest rate r, time to maturity 
r, and volatility a, we have 

E {exp{-rT} {St - K)^) = E (^BS (^So,K,r,T,dr^^ , (5.10) 

that is, we firstly simulate Vsds and ^/VsdWg, using Algorithm [H and 
then compute the Black-Scholes price, for the particular values of So and a cor- 
responding to the trajectory of Vgds and y/V^dW^. This can of course be 
expected to reduce the variance, essentially, we do not estimate the Black-Scholes 
price using Monte Carlo simulation, which is done in Step 4 of Algorithm [H but 
compute the value exactly. Finally, it can be expected that when combining the 
conditional Monte Carlo approach with quasi-Monte Carlo points, the approach 
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Number of simulation trials 


Price estimates 


Standard error estimates 


960 


0.441047 


0.00286267 


1920 


0.443279 


0.00195695 


3840 


0.442995 


0.000698016 


7680 


0.44392 


0.000461543 



Table 2: Price estimates and standard error estimates for a European call option 
using quasi-Monte Carlo points 

is even more efficient, see [29]. This is due to the fact that taking the condi- 
tional expectation has a smoothing effect, which can be expected to improve the 
performance of quasi-Monte Carlo methods, see [22], Subsection 10.1. 

Regarding the quasi-Monte Carlo point sets, we employ the two-dimensional 
Sobol sequence, which is well-known to have the optimal quality parameter t = 0. 
To be more precise, we use two-dimensional Sobol nets, comprised of 2*" points, 
where m = 5, 6, 7, and 8, using the first coordinate to generate Xt and the sec- 
ond coordinate to generate the conditional integrated variance. Furthermore, we 
randomize the nets using Owen's scrambling algorithm, [21], implemented using 
the algorithm presented in [L7j: we produce 30 independent copies of each net, 
allowing us to estimate standard errors. Finally, we remark that besides the abil- 
ity to estimate standard errors, the randomized point sets can also be expected 
to produce better convergence rates, see [11], [25] . 

We use the set of parameters (14. 9 p and, in Table [21 we report estimates of the 
option price and standard errors. The number of simulation trials performed is 
30 * 2™, where m = 5, 6, 7, and 8, as we subject each Sobol net to 30 randomiza- 
tions. 

We note that conditional Monte Carlo combined with quasi-Monte Carlo point 
sets provides us with a substantial variance reduction. Exploring similar variance 
reduction techniques applicable to other payoff functions poses interesting and 
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important future research questions. 



6 Conclusion 



In the present paper, we provided an exact simulation algorithm for the 3/2 
model. A result by Craddock and Lennox allowed us to adapt the Broadie- 
Kaya algorithm for afiine processes to the 3/2 model. Furthermore, we discussed 
variance reduction techniques and found that conditional Monte Carlo combined 
with quasi-Monte Carlo point sets provided significant variance reduction. 

In future work, we aim to discuss path-dependent payoffs under the 3/2 model 
and provide effective variance reduction techniques for path-dependent payoffs. 
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